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Preface 


The  Saugus  River  Floodgate  numerical  modeling  study,  as  documented  in 
this  report,  was  performed  for  the  U.S.  Army  Engineer  Division,  New  England 
(NED).  Mr.  C.  J.  Wener,  Chief  of  NED’s  Hydraulics  and  Water  Quality 
Branch,  was  point  of  contact. 

The  study  was  conducted  in  the  Hydraulics  Laboratory  (HL)  of  the 
U.S.  Army  Engineer  Waterways  Experiment  Station  (WES)  during  the  period 
September  1990  to  June  1992  under  the  direction  of  Messrs.  F.  A.  Herrmann, 
Jr.,  Director,  HL;  R.  A.  Sager,  Assistant  Director,  HL;  W.  H.  McAnally,  Jr., 
Chief,  Estuaries  Division  (ED),  HL;  and  W.  D.  Martin,  Chief,  Estuarine 
Engineering  Branch  (EEB),  ED. 

This  study  was  conducted  by  Dr.  Hsin-Chi  J.  Lin,  EEB,  and  the  report  was 
prepared  by  Dr.  Lin  and  Mr.  David  R.  Richards,  Chief,  Estuarine  Simulation 
Branch  (ESB),  ED. 

Dr.  Robert  W.  Whalin  was  Director  of  WES  during  the  publication  of  this 
report.  COL  Leonard  G.  Hassell,  EN,  was  Commander. 


Conversion  Factors, 
Non-SI  to  SI  Units  of 
Measurement 


Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI  units 
as  follows: 
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1  Introduction 


Background 

The  Saugus  and  Pines  River  estuary  is  located  along  the  Atlantic  coast 
approximately  10  miles1  north  of  Boston,  MA,  near  the  cities  of  Lynn, 

Malden,  and  Revere,  and  the  town  of  Saugus  (Figure  1).  The  Saugus  and 
Pines  Rivers  and  their  tributaries  compose  a  47-square-mile  watershed  area  that 
drains  into  a  tidal  estuary  at  the  mouths  of  the  rivers.  These  estuaries  and  the 
adjacent  saltwater  marshes  total  approximately  1,660  acres.  Freshwater  flows 
from  the  Saugus  and  Pines  Rivers  are  relatively  small.  Storm  water  drainage 
is  temporarily  stored  in  many  areas  in  the  form  of  surface  ponding  when  the 
tide  is  high  followed  by  drainage  when  the  tide  is  falling. 

Because  of  the  topography  and  hydraulics  of  the  Saugus  and  Pines  river 
basins,  storm  events  that  increase  tide  levels  create  a  significant  potential  for 
flooding.  In  1978,  the  eastern  New  England  coastline  was  struck  by  a  storm 
that  created  a  100-year  tidal  flood  event.  The  storm  caused  widespread  and 
record-setting  flooding  in  residential,  commercial,  and  transportation  areas 
along  the  Saugus  and  Pines  Rivers. 

A  plan  was  developed  by  the  U.S.  Army  Engineer  Division,  New  England, 
to  provide  flood  damage  reduction  against  the  Standard  Project  Northeaster 
(SPN)  event  for  nearly  the  entire  project’s  protected  area.  The  principal 
component  of  this  plan  is  construction  of  tidal  floodgates  at  the  mouth  of  the 
Saugus  River.  These  floodgates  will  prevent  tidal  surges  from  entering  the 
river,  thereby  reducing  flooding  within  the  study  area.  The  floodgates  will  be 
constructed  to  maintain  both  safe  navigation  and  natural  tide  levels  and 
flushing  patterns  in  the  estuary  under  normal  conditions.  The  gates  will  be 
closed  only  when  the  projected  tide  levels  are  expected  to  cause  significant 
damage. 

A  field  investigation  by  the  New  England  Division  showed  that  the 
phragmites  reed,  which  indicates  deterioration  in  saltwater  wetlands,  appears  to 
be  expanding  in  the  northwest  comer  of  the  marsh.  Due  to  Federal,  State,  and 


*  A  table  of  factors  for  converting  non-Sl  units  of  measurement  to  SI  units  is  found  on 
page  v. 
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Figure  1.  Vicinity  and  location  maps 


local  interest  in  preserving  and  restoring  this  wetland  by  breaching  the 
abandoned  1-95  embankment  (Figure  1),  a  breaching  plan  was  developed  by 
the  New  England  Division  to  restore  degraded  wetlands  and  increase  tide 
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levels,  resulting  in  increased  flushing  of  nearly  500  acres.  The  1-95  embank¬ 
ment  is  an  abandoned  highway  fill  remaining  from  roadway  construction  activ¬ 
ity  that  was  never  completed.  The  plan  includes  breaching  the  1-95  embank¬ 
ment  at  the  east  branch  of  the  Pines  River  and  widening  the  existing  Pines 
River  opening  in  the  1-95  embankment. 


Objectives 

The  objectives  of  the  numerical  model  study  were  to 

a.  Describe  the  existing  hydrodynamics  of  the  Saugus  and  Pines  Rivers. 

b.  Provide  upstream  and  downstream  boundary  conditions  for  testing  the 
proposed  floodgate  plan  in  a  physical  model  study. 

c.  Determine  the  impacts  caused  by  breaching  of  the  1-95  embankment  at 
the  east  branch  of  Pines  River  and  widened  Pines  River  openings  in  the 
1-95  embankment. 

d.  Evaluate  the  impacts  of  the  floodgate  structure  on  basin  tide  levels, 
circulation  patterns,  storm  surges,  and  sedimentation  and  the  effect  of 
sea  level  rise  on  these  responses. 


Approach 


A  hybrid  modeling  study  approach  (a  combination  of  physical  and  numeri¬ 
cal  models)  was  chosen  to  address  the  numerous  and  complex  concerns  raised 
by  the  construction  and  operation  of  the  proposed  floodgates  An  undistorted, 
l:50-scale  physical  model  representing  a  section  of  the  estuary  was  con¬ 
structed  It  provided  detailed  three-dimensional  information  needed  for  the 
design  of  the  gate  structure.  The  results  of  the  physical  model  provided  input 
to  the  navigation  study. 

Results  from  the  physical  model  and  navigation  studies  are  documented  by 
Brogdon  (in  preparation)  and  Park  (in  preparation),  respectively. 

The  approach  chosen  for  the  numerical  modeling  portion  of  the  overall 
study  was  to  use  the  TABS-MD  numerical  modeling  system  (Thomas  and 
McAnally  1991)  to  evaluate  the  hydrodynamics  and  sediment  transport  of  the 
proposed  plans  and  compare  these  to  existing  conditions.  A  finite  element 
mesh  that  covered  the  Saugus  and  Pines  estuary  and  Lynn  Harbor  was  con¬ 
structed  to  include  marshy  areas  in  the  upper  Pines  River.  The  designed  mesh 
was  sufficiently  refined  to  allow  the  model  to  handle  the  wetting  and  drying 
process  properly  in  the  marshy  areas  and  to  allow  the  observation  of  circula¬ 
tion  patterns  in  the  floodgate  area.  The  numerical  model  was  also  used  to 
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provide  boundary  condition  information  with  which  to  control  the  physical 
model. 

Historically,  the  Saugus  River  near  the  proposed  floodgate  has  not 
experienced  shoaling  problems  and  it  is  not  likely  that  the  project  will  cause 
new  sedimentation  problems.  However,  to  quantify  this  assessment,  a 
sensitivity  study  using  a  sediment  transport  model  was  performed  to  evaluate 
shoaling  potential. 
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2  Description  of  the  Model 


TABS-MD 


TABS-MD  (Thomas  and  McAnally  1991)  is  the  name  of  a  family  of 
computer  programs  used  in  the  multidimensional  modeling  of  hydrodynamics 
(RMA-2V),  sedimentation  (STUDH),  and  constituent  transport  (RMA-4)  in 
river*,  reservoirs,  bays,  and  estuaries.  In  this  study  only  two-dimensional 
models  were  used.  The  system  contains  all  of  the  necessary  preprocessing  and 
postprocessing  utilities  to  allow  relatively  user-friendly  applications.  A  more 
detailed  description  of  the  models  in  the  TABS-MD  system  appears  in 
Appendix  A. 


Mesh  Design 

A  numerical  model  mesh  was  designed  to  allow  replication  of  tidal 
circulation  throughout  the  Saugus  ai;d  Pines  estuary,  Lynn  Harbor,  and  the 
marshy  areas  in  the  upper  Pines  River  'Figure  2).  Tne  mesh  was  sufficiently 
refined  to  model  the  wetting  "id  drying  response  in  the  marshy  areas.  The 
mesh  was  generated  with  high  resolution  in  the  channel  to  represent  the  flow 
patterns  in  sufficient  detail  to  allow  accurate  representation  of  currents  near  the 
floodgate  areas.  The  mesh  consists  of  more  than  8,OOU  elements  and  19,000 
nodes,  which  were  necessary  to  accurately  represent  the  bathymetry  of  the 
study  area. 

A  total  of  four  meshes  were  developed  over  the  course  of  the  study.  The 
first  mesh  (existing  condition)  was  developed  for  the  existing  condition  of  the 
Saugus  and  Pines  estuary  that  extended  from  Broad  Sound  to  the  mans!  v  areas 
in  the  upper  Pines  River  (Figure  2).  The  second  mesh  (base  condition)  was 
developed  from  the  first  mesh  oy  adding  dements  to  account  for  a  breach 
section  at  the  east  branch  of  the  Pines  River  (where  creeks  were  previously  cut 
off  by  the  1-95  fill)  and  widening  the  Pines  River  opening  in  1-95  (Figure  3). 
The  third  mesh  (Plan  2C+7)  was  developed  from  the  second  mesh  by  adding 
more  elements  to  the  floodgate  area  and  including  the  layout  of  the  structure 
(Figure  4).  The  fourth  mesh  (Plan  3)  was  developed  from  the  third  mesh  by 
deleting  the  mesh  upstream  of  the  floodgate  to  simulate  the  floodgate  closed 
condition  (Figure  5), 
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Figure  2.  Saugus  River  estuary  numerical  model  mesh,  existing  condition  (mesh  1) 


Model  geometry  for  the  existing  conditions  was  defined  by  the  U.S. 
Geological  Survey  (USGS)  l:25,000-scale  metric  topographical  map,  Lynn, 
Massachusetts,  dated  1985,  and  Boston  North,  Massachusetts,  dated  1979. 
Additional  bathymetry  was  provided  by  the  New  England  Division  in  the  areas 
of  poor  coverage.  These  include  the  marshy  areas  of  Diamond  Creek  and  the 
upper  Pines  and  Saugus  Rivers,  and  the  deeper  waters  of  Broad  Sound,  Lynn 
Harbor,  and  the  proposed  floodgate  area. 
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Figure  3.  Saugus  River  estuary  numerical  model  mesh,  base  condition 
(mesh  2) 


Numerical  Hydrodynamic  Model 

Boundary  conditions 

Boundary  conditions  for  the  study  were  obtained  from  a  field  survey 
conducted  by  the  U.S.  Army  Engineer  Waterways  Experiment  Station  between 
30  October  and  8  December  1990.  Tid<.  nd  velocity  data  were  collected  to 
provide  boundary  conditions  for  the  model  as  well  as  to  use  as  a  verification 
data  set.  The  14-hr  intensive  velocity  data  collection  period  on  3  November 
1990  encompassed  an  entire  tide  cycle  during  a  spring  tide.  The  tidal  range 
during  the  survey  was  13.13  ft  Mostly  clear  skies  existed  at  the  time  of  the 
survey,  and  wind  conditions  ranged  from  a  slight  breeze  to  light  winds  of  4  to 
5  mph.  A  detailed  report  on  the  field  survey  is  given  by  Fagerburg  et  al. 
(1991).  The  tide  gauge  S0.6,  located  in  Broad  Sound,  was  used  as  a 
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Figure  4.  Saugus  River  estuary  numerical  model  mesh,  Plan  2C+7  (mesh  3) 


water-surface  elevation  boundary  condition.  The  drainage  areas  of  both 
Saugus  and  Pines  Rivers  are  very  small,  and  less  than  1.0  cfs  of  inflow  was 
observed  in  the  field  trip.  Therefore,  no  freshwater  inflow  boundary  was 
specified  in  the  Saugus  and  Pines  Rivers. 

Model  parameters 

Model  verification  resulted  in  one  final  set  of  model  parameters 
representing  Manning’s  n  values  and  eddy  viscosities.  The  parameter  values 
were  selected  by  adjustment  within  a  range  of  realistic  values  until  an  optimum 
comparison  of  the  model’s  computed  water  levels  and  currents  to  field 
measurements  was  obtained.  The  following  tabulation  lists  eddy  viscosity 
coefficients  and  Manning’s  values  used  in  this  study: 
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Figure  5.  Saugus  River  estuary  numerical  model  mesh,  floodgate  closed 
condition  (mesh  4) 
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Location 

Eddy  Viscosity,  Ib-ssc/tt2 

Manning's  n 

Broad  Sound 

50 

0.020 

Saugus  River 

50-100 

0.020-0  040 

Pines  River 

50 

0.040-0.080 

Marshy  area 

100 

0.070 

Verification 

The  field  verification  data  were  obtained  from  9  water  level  stations  (Fig¬ 
ure  6)  and  11  current  stations  (Figure  7)  located  in  Saugus  and  Pines  Rivers. 
Water  level  elevation  measurements  at  each  station  were  recorded  using  Micro¬ 
tide  water  level  recorders.  Velocities  at  each  station  were  measured  with  the 
deployment  of  recording  instruments.  A  Gurley  Model  665  vertical-axis,  cup- 
type  impeller  velocity  meter  with  direct  velocity  readout  capabilities  was  used 
to  measure  current  speeds.  At  each  station  the  velocity  data  were  measured  at 
three  depths — surface,  middepth,  and  bottom — for  each  hour  of  the  survey 
period.  The  bottom  measurement  was  made  2  ft  from  the  actual  bottom.  The 
middepth  data  were  obtained  at  the  calculated  middepth.  The  surface  measure¬ 
ment  was  obtained  2  ft  below  the  water  surface. 

Because  of  high  tidal  amplitudes,  RMA-2V  was  operated  with  a  variable 
time-step  ranging  from  5  to  15  min.  The  simulation  started  when  the  tide  was 
at  its  peak  in  Broad  Sound.  A  6-hr  initial  period  allowed  the  transients 
induced  by  initialization  (spin-up)  of  RMA-2V  to  dissipate  and  the  model 
solution  to  respond  correctly  to  imposed  boundary  conditions.  A  1-day  data 
set  (two  tidal  cycles)  was  used  for  tide  verification.  A  14-hr  period  of  velocity 
data  was  used  for  velocity  verification. 


Tide  verification  results 

Water  levels  from  the  mode!  and  field  measurements  are  compared  in  Plates 
1-8.  The  water  level  comparisons  indicated  good  agreement  between  model 
results  and  field  measurements.  The  entire  estuary  is  nearly  in  phase  in  both 
model  and  field  data  except  at  tide  Gages  S4.4,  S9.1,  and  S9.3.  All  three 
stations  are  located  in  marshy  areas.  At  Gage  S4.4,  located  in  Diamond  Creek, 
comparisons  of  rising  and  falling  tide  calculated  by  the  model  with  the  field 
measurements  showed  good  agreement;  however,  the  computed  peak  water 
level  was  slightly  high  (about  0.3  ft  at  peak  tide)  and  the  phasing  of  peak  tide 
was  off  by  about  15  min.  At  Gages  S9.1  and  S9.3,  located  in  the  upper  Pines 
River,  the  phasing  of  flood  tide  was  good;  however,  the  computed  water  levels 
were  slightly  high  (about  0.2  to  0.3  ft)  at  the  peak  tide  and  the  ebb  tide  fell 
faster  than  the  field  data  by  about  1.0  hr.  The  marshes  are  full  of  dense 
grasses  (about  1  ft  high)  and  many  interconnected  mosquito  ditches  (1  to  2  ft 
wide  and  2  to  4  ft  deep).  Water  was  stored  in  the  marshes  and  ditches  during 
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Figure  6.  Water  level  field  stations 


the  high  tide  and  was  drained  slowly  to  the  channel  during  the  low  tide.  At 
Gage  S9.5,  located  in  upper  Saugus  River,  the  computed  water  levels  were 
slightly  high  (about  0.2  ft)  at  the  peak  ebb  tide  and  the  phasing  was  off  by 
about  1  hr.  In  general,  the  computed  water  levels  showed  good  agreement, 
except  in  the  marshy  areas.  The  computed  water  level  in  the  marshy  area  was 
slightly  high  by  about  0.3  ft.  This  level  of  verification  is  the  best  possible 
given  the  present  state  of  technology.  A  more  highly  refined  mesh  that 
included  each  ditch  would  have  provided  better  verification,  but  it  was  not 
possible  to  include  all  ditches.  The  tidal  verification  is  acceptable  for  the 
purpose  of  this  study. 
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Figure  7.  Ranges  and  locations  for  14-hr  velocity  data  collection 


Velocity  verification  results 

Velocities  predicted  by  the  model  are  compared  with  field-measured 
values  in  Plates  9-19.  At  all  survey  stations,  the  vertical  depth  profile  from  the 
field  survey  was  averaged  to  obtained  a  single  value  for  comparison  with  the 
depth-averaged  model.  The  model-predicted  velocities  compared  quite 
favorably  with  the  vertically  averaged  field  measurements.  Small  variations  of 
model  results  compared  to  field  measurements  occurred  at  sta  1.0A  and  1.0B 
in  Broad  Sound.  Good  comparisons  occurred  at  sta  2.0A  and  2.0B,  located 
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near  the  proposed  floodgates  area.  At  sta  3.0B,  4.0B,  5.0B,  and  6.0B,  located 
in  the  Pines  River,  peak  flood  velocities  were  underpredicted  by  about  0.4  fps 
and  peak  ebb  velocities  were  overpredicted  by  about  1.0  fps.  Good 
comparisons  occurred  at  all  stations  in  the  Saugus  River,  except  7.0B,  where 
peak  flood  velocities  were  underpredicted  by  about  0.4  fps  and  peak  ebb 
velocities  were  overpredicted  by  about  0.4  fps.  All  of  the  underprediction  in 
peak  flood  velocities  and  overprediction  in  peak  ebb  velocities  can  be 
explained  by  model  elements  being  removed  or  added  when  an  element 
became  dry  or  wet.  When  the  water  level  falls  and  the  element  becomes  dry 
or  when  the  water  level  rises  and  the  element  becomes  wet  again,  this 
condition  causes  numerical  instability  and  shocks  the  system.  The  numerical 
shocks  caused  by  wetting  and  drying  are  one  of  the  most  difficult  problems  in 
numerical  modeling.  It  can  be  improved  by  adding  more  resolution  in  the 
areas  when  the  wetting  and  drying  occur,  but  the  penalty  is  more  computer 
time.  The  level  of  accuracy  in  the  velocity  verification  was  as  good  as  could 
be  obtained  within  time  and  cost  constraints  and  is  considered  adequate  for  the 
purposes  of  the  study. 

Plates  20  and  21  provide  model  velocity  vector  plots  for  maximum  flood 
and  maximum  ebb,  respectively,  for  the  existing  conditions.  These  plots  show 
the  general  flow  pattern  in  the  Saugus  and  Pines  Rivers. 


Numerical  Sediment  Transport  Model 

Boundary  conditions  for  the  STUDH  sediment  transport  model  consisted  of 
suspended  sediment  concentrations.  The  nodal  velocities  from  the  hydro- 
dynamic  model  were  saved  and  used  to  update  the  velocity  field  in  the  model 
at  the  beginning  of  each  time-step. 

Since  the  floodgate  area  has  not  experienced  sediment  problems,  the  sedi¬ 
ment  study  focused  on  a  sensitivity  analysis  of  model  parameters.  The 
STUDH  code  uses  the  Ackers-White  sediment  transport  function  for  computing 
noncohesive  sediment  transport  (Ackers  and  White  1973). 

Sediment  data  were  insufficient  to  verify  the  sediment  model.  However, 
some  suspended  sediment  concentration  data  samples  were  collected  at  the 
individual  sampling  stations  during  the  14-br  survey  (Fagerburg  et  al.  1991). 
The  suspended  sediment  concentrations  were  found  to  be  low  during  the  sur¬ 
vey  period.  The  suspended  sediment  concentrations  at  sta  R1.0A  were  less 
than  5  mg /C  during  the  flood  tide  and  about  30.0  mg ft  during  the  ebb  tide.  It 
appears  that  the  sediment  concentration  in  the  study  areas  is  very  low  during 
normal  tide  conditions.  A  constant  suspended  sediment  concentration  of 
30  mg/C  was  specified  at  the  boundary  in  order  to  produce  some  shoaling  for 
an  evaluation  of  parameter  sensitivity.  The  various  input  variables  included 
diffusion  coefficients,  Manning’s  n  to  compute  bed  shear  stress,  and  effective 
particle  size  for  transport.  Each  was  adjusted  until  the  model  produced  reason¬ 
able  shoaling  and  scour  patterns.  The  model  parameter  values  are  listed  in  the 
following  tabulation. 
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Parameter 

Value 

Effective  particle  size  for  transport,  mm 

0.10 

Effective  setting  velocity,  m/sec 

0.007 

Diffusion  coefficient.  m2/sec 

10 

Manning's  n  for  bed  shear  stress 

0.045 

Computation  time-step,  min 

15 

Plate  22  shows  the  deposition  and  erosion  pattern  after  a  24-hr  simulation 
near  the  proposed  floodgate  area.  The  simulation  indicated  that  there  is  little 
deposition  or  erosion  near  the  proposed  floodgate  area.  The  deposition 
occurred  at  the  boundary  in  Broad  Sound  and  was  caused  by  the  excessive 
sediment  load  (30  mg/i)  at  the  boundary. 

The  sensitivity  analysis  of  sediment  movement  was  conducted  for  the 
normal  tide  conditions.  It  did  not  address  the  sediment  movement  associated 
with  runoff,  storm  surge,  and  wave-producing  storm  events.  The  analysis  was 
focused  on  any  change  of  the  sediment  deposition  and  scour  pattern  under  the 
proposed  floodgate  as  compared  to  the  existing  condition.  A  24-hr  simulation 
(about  two  tidal  cycles)  was  used  to  indicate  any  significant  change  in 
sediment  deposition  and  scour  pattern  in  the  estuary. 
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3  Model  Results 


Hydrodynamics 

The  verified  hydrodynamic  model  was  used  to  provide  water  levels  and 
discharges  as  the  boundary  conditions  for  use  in  the  physical  model.  It  was 
also  used  to  evaluate  the  impacts  of  breaching  the  1-95  embankment  and  to 
determine  the  effects  of  the  proposed  floodgate,  rising  sea  levels,  and  storm 
surges  on  the  area. 


Physical  model  boundary  conditions 

The  computed  water  levels  and  discharges  from  the  numerical  model  at 
cross  sections  1,  2,  3,  4,  and  5  (Figure  8)  were  provided  for  use  in  the  physical 
model  study.  Two  sets  of  water  levels  and  discharges  for  existing  and  base 
conditions  were  based  on  the  maximum  discharges  (both  flood  and  ebb  tide)  at 
the  proposed  floodgate  site  during  spring  and  neap  tide  cycles. 

Table  1  shows  the  boundary  conditions  for  the  physical  model  under  the 
existing  spring  tide  conditions.  The  flows  were  29,000  and  22,000  cfs  for  the 
flood  and  ebb  tides,  respectively.  The  water  levels  were  3.48  ft1  and  4.81  ft 
for  the  flood  and  ebb  tides,  respectively. 

Table  2  shows  the  boundary  conditions  for  the  physical  model  under  the 
existing  neap  tide  conditions.  The  flows  were  12,300  and  11,300  cfs  for  the 
flood  and  ebb  tides,  respectively.  The  water  levels  were  1.46  ft  and  -0.30  ft 
for  the  flood  and  ebb  tides,  respectively. 

Table  3  shows  the  boundary  conditions  for  the  physical  model  under  the 
base  spring  tide  conditions.  The  flows  were  27,100  and  28,200  cfs  for  the 
flood  and  ebb  tides,  respectively.  The  water  levels  were  6.24  ft  and  2.70  ft  for 
the  flood  and  ebb  tides,  respectively. 


*  All  elevations  (el)  cited  herein  are  in  feet  referred  to  the  National  Geodetic  Vertical  Datum 
(NGVD). 
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Figure  8.  Locations  of  boundary  conditions  for  physical  model  study 


Table  4  shows  the  boundary  conditions  for  the  physical  model  under  the 
base  neap  tide  conditions.  The  flows  were  1 1,900  and  10,900  cfs  for  the  flood 
and  ebb  tides,  respectively.  The  water  levels  were  1.46  ft  and  -0.31  ft  for  the 
flood  and  ebb  tides,  respectively. 

Breaching  of  1*95  embankment 

The  purpose  of  considering  breaching  the  abandoned  1-95  embankment  and 
widening  the  Pines  River  opening  through  the  embankment  was  to  restore 
deteriorated  wetlands  in  the  Pines  River.  If  the  breaching  was  implemented,  it 
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was  important  to  determine  the  increase  in  tide  level  in  the  upper  Pines  River 
due  to  the  cut  of  1-95.  Tlie  layout  and  dimensions  of  the  breach  site 
(700  sq  ft)  in  1-95  and  the  enlarged  opening  (4,200  sq  ft)  in  Pines  River 
(Figure  9)  were  provided  by  the  New  England  Division. 

After  the  abandoned  1-95  embankment  was  breached,  the  Route  107  bridge 
opening  became  a  control  section.  The  existing  opening  is  very  small  and  the 
computed  velocity  was  more  than  6.0  fps  (Plate  23).  A  2-fps  target  velocity 
was  suggested  by  New  England  Division.  Two  flow  areas  (40  ft  wide  and 
90  ft  wide)  at  the  Route  107  bridge  were  studied  based  on  the  2-fps  criteria. 
Plates  24  and  25  show  that  the  velocity  at  the  bridge  will  be  about  3.0  fps  for 
a  40-ft-wide  opening  and  1.5  fps  for  a  90-ft-wide  opening,  respectively. 

Plates  26-28  show  that  the  water  levels  in  the  marshes  will  not  change 
significantly  (less  than  0.05  ft). 

Plates  29-36  compare  the  computed  water  levels  for  the  existing  and  base 
conditions.  The  breaching  of  the  1-95  embankment  will  increase  the  water 
level  in  marshy  areas  of  the  upper  Pines  River  by  about  0.5  ft  at  the  peak  tide 
during  the  spring  tide  condition  and  will  not  affect  the  water  level  in  the 
marshy  area  of  Diamond  Creek.  The  time  lag  between  the  peak  water  levels  at 
sta  S0.6  and  at  sta  S9.3  is  about  2  hr  in  the  existing  condition.  The  time  lag 
was  reduced  to  1.0  hr  due  to  breaching  of  the  abandoned  1-95  embankment. 


Floodgate 

The  layout  and  dimensions  of  the  proposed  floodgate  (Plan  2C+7)  and 
approach  channel  were  provided  by  the  physical  model  study  (Figure  10).  The 
third  mesh  (Figure  4)  with  the  same  boundary  conditions  as  specified  for  exist¬ 
ing  conditions  was  the  input  to  the  model.  Plates  37-44  compare  the  computed 
water  levels  for  the  existing  condition  and  the  proposed  floodgate.  The  results 
indicated  no  measurable  change  in  water  levels  in  the  estuary  except  in  the 
marshy  area  of  the  upper  Pines  River.  Water  levels  at  sta  S9.1  and  S9.3  will 
increase  about  0.5  ft  at  the  peak  tide  during  the  spring  tide  condition,  and  the 
time  to  peak  water  level  will  decrease  about  1.0  hr  compared  to  the  existing 
condition.  Plates  45-52  compare  the  computed  water  levels  for  the  base  condi¬ 
tion  and  the  proposed  floodgate.  The  results  show  no  measurable  change  in 
estuary  water  levels  due  to  the  floodgate  compared  to  the  base  condition. 

Based  on  these  results,  it  is  concluded  that  floodgates  will  not  cause  significant 
change  of  water  levels  in  the  Pines  and  Saugus  Rivers.  Plates  53-55  compare 
the  computed  velocity  for  the  base  condition  and  the  Plan  2C+7  condition  at 
the  enlarged  opening  in  the  Pines  River,  at  the  breach  cross  section,  and  at  the 
upstream  side  of  the  General  Edwards  bridge  in  Saugus  River.  The  results 
show  no  significant  change  of  velocity.  Plates  56  and  57  show  the  model 
velocity  vector  plots  for  maximum  flood  and  ebb,  respectively.  Plates  58  and 
59  show  the  model  velocity  vector  plots  for  maximum  flood  and  ebb,  respec¬ 
tively,  near  the  floodgate  area. 
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Sea  level  rise 


The  third  mesh  was  used  for  this  task  with  1.0-fl  elevated  tide  boundary 
conditions.  Plates  60-67  compare  the  computed  water  levels  between  the 
spring  tide  and  elevated  spring  tide  under  Plan  2C+7.  The  water  levels  at  the 
peak  flood  tide  and  the  peak  ebb  tide  will  increase  1.0  ft  in  the  study  area. 

The  effects  of  water  level  changes  due  to  the  1  -ft  rise  in  sea  level  for  the  exist¬ 
ing  and  the  base  condition  were  not  included  in  the  study.  Based  on  the 
results  discussed  in  the  sections  "Breaching  of  1-95  embankment"  and  "Rood- 
gate,"  the  floodgate  will  not  significantly  affect  the  water  levels  in  the  estuary 
for  the  1-ft  rise  in  sea  level. 


Storm  surges 

The  floodgate  will  be  closed  when  the  projected  tide  levels  are  expected  to 
cause  significant  damage.  The  fourth  mesh  with  the  same  boundary  conditions 
as  specified  in  existing  conditions  was  the  input  to  the  model.  Comparisons  of 
the  computed  water  levels  ir.  the  proposed  floodgate  area  between  the  existing 
and  the  floodgate  closed  con.'  tton  show  no  measurable  differences  (Plate  68). 
The  results  indicated  that  the  closure  of  the  floodgate  will  not  cause  significant 
change  of  water  level  in  Broad  Sound.  This  comparison  is  based  on  hydro¬ 
dynamics  and  does  not  include  effects  of  wind  and  wave  setup.  This  analysis 
was  just  a  sensitivity  test  to  see  if  closure  of  the  floodgate  would  cause  any 
change  in  tide  level  in  Broad  Sound. 


Sedimentation 

Plate  69  shows  the  shoaling  and  scour  pattern  near  the  proposed  floodgate 
area  under  Plan  2C+7.  The  plot  was  generated  based  on  the  same  input  data 
as  specified  in  the  existing  condition  (Plate  22).  The  results  indicated  little 
difference  in  deposition  and  erosion  pattern  compared  to  the  existing  condition. 
The  proposed  floodgate  will  not  alter  the  shoaling  and  scour  pattern  in  the 
study  area,  but  local  shoaling  and  scour  may  occur  near  the  proposed  floodgate 
pier. 


20 


Chapter  3  Model  Results 


4  Conclusions 


The  RMA-2V  model  was  successfully  verified  to  limited  field  measure¬ 
ments  including  a  14-hr  field  survey  of  water  levels  and  velocity  measure¬ 
ments.  The  comparisons  of  the  computed  water  levels  and  velocities  to  field 
measurements  were  good.  At  the  stations  in  marshy  areas,  the  computed  water 
levels  were  slightly  high  by  about  0.2  ft  and  the  ebb  tide  fell  about  1  hr  faster 
than  the  field  data.  The  model  error  was  suspected  to  be  the  result  of  difficul¬ 
ties  in  properly  representing  the  storage  of  water  in  marshes  and  interconnected 
mosquito  ditches.  In  the  Pines  River,  the  velocities  were  underpredicted  by 
about  0.3  fps  for  peak  flood  and  overpredicted  about  1.0  fps  for  peak  ebb. 
These  variations  were  caused  by  the  elements  that  are  removed  when  they  are 
dry  and  added  when  they  are  wet  again.  These  differences  were  small  and  not 
expected  to  significantly  impact  the  use  of  model  results  for  the  intended 
purposes. 

The  study  provided  boundary  conditions  for  the  physical  model  study 
under  the  existing  and  base  conditions  for  both  spring  and  neap  tides.  For  the 
existing  spring  tide  conditions,  the  flows  were  29,000  and  22,000  cfs  for  the 
flood  and  ebb  tides,  respectively.  The  corresponding  water  levels  were  3.48  ft 
and  4.81  ft  for  the  flood  and  ebb  tides,  respectively.  For  the  existing  neap  tide 
conditions,  the  flows  were  12,300  and  11,300  cfs  for  the  flood  and  ebb  tides, 
respectively.  The  corresponding  water  levels  were  1.46  ft  and  -0.30  ft  for  the 
flood  and  ebb  tides,  respectively.  For  the  base  spring  tide  conditions,  the 
flows  were  27,100  and  28,200  cfs  for  the  flood  and  ebb  tides,  respectively. 

The  corresponding  water  levels  were  6.24  ft  and  2.70  ft  for  the  flood  and  ebb 
tides,  respectively.  For  the  base  neap  tide  conditions,  the  flows  were  11,900 
and  10.900  cfs  for  the  flood  and  ebb  tides,  respectively.  The  corresponding 
water  levels  were  1.46  ft  and  -0.31  ft  for  the  flood  and  ebb  tides,  respectively. 

Breaching  of  the  abandoned  1-95  embankment  and  widening  the  Pines 
River  opening  on  1-95  will  increase  tidal  flow  in  marshy  areas.  The  water 
levels  in  marshy  areas  will  increase  about  0.5  ft  at  the  peak  tide  under  a  spring 
tide  condition.  The  time  lag  of  the  peak  water  levels  between  the  Broad 
Sound  and  upper  marshy  areas  was  reduced  from  2  hr  to  1  hr. 

Plan  2C+7  will  not  cause  significant  change  of  water  levels  in  the  Pines 
and  Saugus  Rivers  under  the  normal  tide  conditions.  It  will  protect  the  study 
areas  from  flooding  during  the  storm  events. 
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The  water  levels  in  the  marshy  areas  under  Plan  2C+7  will  increase  about 
1.0  ft  at  the  peak  flood  tide  and  ebb  tide  for  the  1-ft  rise  in  sea  level. 

The  proposed  floodgate  will  not  alter  the  sediment  deposition  or  scour 
pattern  in  the  estuary  under  the  normal  tide  condition,  but  local  scour  near  the 
piers  may  occur. 
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Total  28,240 


Table  2 

Boundary  Conditions  for  Physical  Model  Under  Existing 

Neap  Tide 

Boundary  Saction 

Water  Surface  Elevation 

Dlacharge,  cfa 

Flood  Direction 

Inflow 

3 

1.462 

5,770 

4 

1.466 

4.200 

5 

1.467 

2,350 

Total 

12,320 

Outflow 


1 

1.412 

5,976 

2 

1.416 

2,995 

Total 

8,991 

Ebb  Direction 

Inflow 

1 

-0.252 

5,275 

2 

-0.248 

2,730 

Total 

8,005 

Outflow  ] 

3 

-0.313 

5,880 

4 

-0.309 

4,160 

5 

•0.306 

1,270 

Total  11.310 


Table  3 

Boundary  Conditions  for  Physical  Model  Under  the  Base 
Spring  Tide 

Boundary  Section  Water  Surface  Elevation  Discharge,  cfa 


2  6.18  7.600 

Total  22,100 


Ebb  Direction 


Inflow 

1  2.92  17,400 


Total  28,200 


Table  4 

Boundary  Conditions  for  Physical  Model  Under  Base 

Neap  Tide 

Boundary  Saction 

Water  Surfaca  Elevation 

Dlecharge,  eta 

Flood  Dirac tion 

Inflow 

3 

1.462 

5,660 

4 

1.466 

3,990 

5 

1.467 

2,220 

Total 

11,870 

j  Outflow  | 

i 

1.419 

5,350 

2 

1.423 

2.960 

Total 

8,330 

Ebb  Direction 

Inflow 


1 

•0.253 

4.940 

2 

-0.251 

2,690 

Total  7,630 

Outflow 

3 

-0.313 

5,740 

4 

-0.309 

3,990 

5 

-0.306 

1,210 

Total  10,940 
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Appendix  A 

The  TABS-MD  System 


TABS-MD  is  a  collection  of  generalized  computer  programs  and  utility 
codes  integrated  into  a  numerical  modeling  system  for  studying 
two-dimensional  hydrodynamics,  sedimentation,  and  transport  problems  in 
rivers,  reservoirs,  bays,  and  estuaries.  A  schematic  representation  of  the  sys¬ 
tem  is  shown  in  Figure  Al.  It  can  be  used  either  as  a  stand-alone  solution 
technique  or  as  a  step  in  the  hybrid  modeling  approach.  The  basic  concept  is 
to  calculate  water-surface  elevations,  current  patterns,  sediment  erosion,  trans¬ 
port  and  deposition,  the  resulting  bed  surface  elevations,  and  the  feedback  to 
hydraulics.  Existing  and  proposed  geometry  can  be  analyzed  to  determine  the 
impact  on  sedimentation  of  project  designs  and  to  determine  the  impact  of 
project  designs  on  salinity  and  on  the  stream  system.  The  system  is  described 
in  detail  by  Thomas  and  McAnally  (1985). 

The  three  basic  components  of  the  system  are  as  follows: 

a.  "A  Two-Dimensional  Model  for  Free  Surface  Flows,"  RMA-2V. 

b.  "Sediment  Transport  in  Unsteady  2-Dimensional  Flows,  Horizontal 
Plane,"  STUDH. 

c.  "Two-Dimensional  Finite  Element  Program  for  Water  Quality,"  RMA-4. 

RMA-2V  is  a  finite  element  solution  of  the  Reynolds  form  of  the 
Navier-Stokes  equations  for  turbulent  flows.  Friction  is  calculated  with 
Manning’s  equation  and  eddy  viscosity  coefficients  are  used  to  define  the 
turbulent  losses.  A  velocity  form  of  the  basic  equation  is  used  with  side 
boundaries  treated  as  either  slip  or  static.  The  model  automatically  recognizes 
dry  elements  and  corrects  the  mesh  accordingly.  Boundary  conditions  may  be 
water-surface  elevations,  velocities,  or  discharges  and  may  occur  inside  the 
mesh  as  well  as  along  the  edges. 

The  sedimentation  model,  STUDH,  solves  the  convection-diffusion  equation 
with  bed  source  terms.  These  terms  are  structured  for  either  sand  or  cohesive 
sediments.  The  Ackers- White  (1973)  procedure  is  used  to  calculate  a  sediment 
transport  potential  for  the  sands  from  which  the  actual  transport  is  calculated 
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Figure  A1 .  TABS-2  schematic 

based  on  availability.  Clay  erosion  is  based  on  work  by  Partheniades  (1962) 
and  Ariathurai  and  the  deposition  of  clay  utilizes  Krone’s  equations 
(Ariathurai,  MacArthur,  and  Krone  1977).  Deposited  material  forms  layers,  as 
shown  in  Figure  A2,  and  bookkeeping  allows  up  to  10  layers  at  each  node  for 
maintaining  separate  material  types,  deposit  thickness,  and  age.  The  code  uses 
the  same  mesh  as  RMA-2V. 

Salinity  calculations,  RMA-4,  are  made  with  a  form  of  the  convective- 
diffusion  equation  which  has  general  source-sink  terms.  Up  to  seven  conserva¬ 
tive  substances  or  substances  requiring  a  decay  term  can  be  routed.  The  code 
uses  the  same  mesh  as  RMA-2V. 

Each  of  these  generalized  computer  codes  can  be  used  as  a  stand-alone 
program,  but  to  facilitate  the  preparation  of  input  data  and  to  aid  in  analyzing 
results,  a  family  of  utility  programs  was  developed  for  the  following  purposes: 

a.  Digitizing 

b.  Mesh  generation 

c.  Spatial  data  management 

d.  Graphical  output 

e.  Output  analysis 

/  File  management 

g.  Interfaces 

h.  Job  control  language 
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Figure  A 2.  Two-dimensional  finite  element  mesh 


Finite  Element  Modeling 

The  TABS-2  numerical  models  used  in  this  effort  employ  the  finite  element 
method  to  solve  the  governing  equations.  To  help  those  who  are  unfamiliar 
with  the  method  to  better  understand  this  report,  a  brief  description  of  the 
method  is  given  here. 
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The  finite  element  method  approximates  a  solution  to  equations  by  dividing 
the  area  of  interest  into  smaller  subareas,  which  are  called  dements  The 
dependent  variables  (e.g.,  water-surface  elevations  and  sediment  concentra¬ 
tions)  are  approximated  over  each  element  by  continuous  functions  which 
interpolate  in  terms  of  unknown  point  (node)  values  of  the  variables.  An  error, 
defined  as  the  deviation  of  the  approximation  solution  from  the  correct  solu¬ 
tion,  is  minimized.  Then,  when  boundary  conditions  arc  imposed,  a  set  of 
solvable  simultaneous  equations  is  created.  The  solution  is  continuous  over 
the  area  of  interest. 

In  one-dimensional  problems,  elements  are  line  segments.  In  twodimen¬ 
sional  problems,  the  elements  are  polygons,  usually  either  triangles  or  quad¬ 
rilaterals.  Nodes  are  located  on  the  edges  of  elements  and  occasionally  inside 
the  elements.  The  interpolating  functions  may  be  linear  or  higher  order 
polynomials.  Figure  A 2  illustrates  a  quadrilateral  element  with  eight  nodes 
and  a  linear  solution  surface  where  F  is  the  interpolating  function. 

Most  water  resource  applications  of  the  finite  element  method  use  the 
Galcrkin  method  of  weighted  residuals  to  minimize  error.  In  this  method  the 
residual,  the  total  error  between  the  approximate  and  correct  solutions,  is 
weighted  by  a  function  that  is  identical  with  the  interpolating  function  and  then 
minimized.  Minimization  results  in  a  set  of  simultaneous  equations  in  terms  of 
nodal  values  of  the  dependent  variable  (e.g.  water-surface  elevations  or  sedi¬ 
ment  concentration).  The  time  portion  of  time-dependent  problems  can  be 
solved  by  the  finite  clement  method,  but  it  is  generally  more  efficient  to 
express  derivatives  with  respect  to  lime  in  finite  difference  form. 


The  Hydrodynamic  Model,  RMA-2V 

Applications 

This  program  is  designed  for  lar-field  problems  in  which  vertical  accelera¬ 
tions  arc  negligible  and  the  velocity  vectors  at  a  node  generally  point  in  the 
same  directions  over  the  entire  depth  of  the  water  column  at  any  instant  of 
time.  It  expects  a  homogeneous  fluid  with  a  free  surface.  Both  steady  and 
unsteady  state  problems  can  be  analyzed.  A  surface  wind  stress  can  be 
imposed. 

The  program  has  been  applied  to  calculate  flow  distribution  around  islands; 
flow  at  bridges  having  one  or  more  relief  openings,  in  contracting  and  expand¬ 
ing  reaches,  into  and  out  of  off-channel  hydropower  plants,  at  river  junctions, 
and  into  and  out  of  pumping  plant  channels;  and  general  flow’  patterns  in 
rivers,  reservoirs,  and  estuaries. 
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Limitations 


This  program  is  not  designed  for  near-field  problems  where  flowstructure 
interactions  (such  as  vortices,  vibrations,  or  vertical  accelerations)  are  of 
interest.  Areas  of  vertically  stratified  flow  are  beyond  this  program’s  cap¬ 
ability  unless  it  is  used  in  a  hybrid  modeling  approach.  It  is  two-dimensional 
in  the  horizontal  plane,  and  zones  where  the  bottom  current  is  in  a  different 
direction  from  the  surface  current  must  be  analyzed  with  considerable  subjec¬ 
tive  judgment  regarding  long-term  energy  considerations.  It  is  a  free-surface 
calculation  for  subcritical  flow  problems. 


Governing  equations 

The  generalized  computer  program  RMA-2V  solves  the  depth-integrated 
equations  of  fluid  mass  and  momentum  conservation  in  two  horizontal  direc¬ 
tions.  The  form  of  the  solved  equations  is 
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where 

h  -  depth 

u,v  =  velocities  in  the  Cartesian  directions 
x,y,t  -  Cartesian  coordinates  and  time 
p  =  density 
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e  =  eddy  viscosity  coefficient,  for  xx  =  normal  direction  on  x-axis 
surface;  yy  -  normal  direction  on  y-axis  surface;  xy  and 
yx  =  shear  direction  on  each  surface 
g  =  acceleration  due  to  gravity 
a  =  elevation  of  bottom 
n  =  Manning’s  n  value 
1.486  =  conversion  from  SI  to  non-Sl  units 
£  =  empirical  wind  shear  coefficient 
Va  =  wind  speed 
Tp  =  wind  direction 
u>  =  rate  of  earth’s  angular  rotation 
<j>  =  local  latitude 

Equations  Al,  A2,  and  A3  are  solved  by  the  finite  element  method  using 
Galerkin  weighted  residuals.  The  elements  may  be  either  quadrilaterals  or 
triangles  and  may  have  curved  (parabolic)  sides.  The  shape  functions  are 
quadratic  for  flow  and  linear  for  depth.  Integration  in  space  is  performed  by 
Gaussian  integration.  Derivatives  in  time  are  replaced  by  a  nonlinear  finite 
difference  approximation.  Variables  are  assumed  to  vary  over  each  time  inter 
val  in  the  form 


/(*)  =  /(0)  ♦  at  *  btc  tQ  s  t  <  t 


(A4) 


which  is  differentiated  with  respect  to  time,  and  cast  in  finite  difference  form. 
Letters  a,  b,  a^^  c  are  constants.  It  has  been  found  by  experiment  that  the  best 
value  for  c  is  1.5  (Norton  and  King  1977). 

The  solution  is  fully  implicit  and  the  set  of  simultaneous  equations  is 
solved  by  Newton-Raphson  iteration.  The  computer  code  executes  the  solution 
by  means  of  a  front-type  solver  that  assembles  a  portion  of  the  matrix  and 
solves  it  before  assembling  the  next  portion  of  the  matrix.  The  front  solver’s 
efficiency  is  largely  independent  of  bandwidth  and  thus  does  not  require  as 
much  care  in  formation  of  the  computational  mesh  as  do  traditional  solvers. 

The  code  RMA-2V  is  based  on  the  earlier  version  RMA-2  (Norton  and 
King  1977)  but  differs  from  it  in  several  ways.  It  is  formulated  in  terms  of 
velocity  (v)  instead  of  unit  discharge  (v/i),  which  improves  some  aspects  of  the 
code’s  behavior;  it  permits  drying  and  wetting  of  areas  within  the  grid;  and  it 
permits  specification  of  turbulent  exchange  coefficients  in  directions  other  than 
along  the  x-  and  2-axes.  For  a  more  complete  description,  see  Appendix  F  of 
Thomas  and  McAnally  (1985). 
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The  Sediment  Transport  Model,  STUDH 

Applications 

STUDH  can  be  applied  to  clay  and/or  sand  bed  sediments  where  flow 
velocities  c  n  be  considered  two-dimensional  (i.e.,  the  speed  and  direction  can 
be  satisfactorily  represented  as  a  depth-averaged  velocity).  It  is  useful  for  both 
deposition  and  erosion  studies  and,  to  a  limited  extent,  for  stream  width 
studies.  The  program  treats  two  categories  of  sediment:  noncohesive,  which 
is  referred  to  as  sand  here,  and  cohesive,  which  is  referred  to  as  clay. 


Limitations 

Both  clay  and  sand  may  be  analyzed,  but  the  model  considers  a  single, 
effective  grain  size  for  each  and  treats  each  separately.  Fall  velocity  must  be 
prescribed  along  with  the  water-surface  elevations,  x-velocity,  y-velocity,  diffu¬ 
sion  coefficients,  bed  density,  critical  shear  stresses  for  erosion,  erosion  rate 
constants,  and  critical  shear  stress  for  deposition. 

Many  applications  cannot  use  long  simulation  periods  because  of  their 
computation  cost.  Study  areas  should  be  made  as  small  as  possible  to  avoid 
an  excessive  number  of  elements  when  dynamic  runs  are  contemplated  yet 
must  be  large  enough  to  permit  proper  posing  of  boundary  conditions.  The 
same  computation  time  interval  must  be  satisfactory  for  both  the  transverse  and 
longitudinal  flow  directions. 

The  program  does  not  compute  water-surface  elevations  or  velocities;  there¬ 
fore  these  data  must  be  provided.  For  complicated  geometries,  the  numerical 
model  for  hydrodynamic  computations,  RMA-2V,  is  used. 


Governing  equations 

The  generalized  computer  program  STUDH  solves  the  depth-integrated 
convection-dispersion  equation  in  two  horizontal  dimensions  for  a  single  sedi 
ment  constituent.  For  a  more  complete  description,  see  Appendix  G  of 
Thomas  and  McAnally  (1985).  The  form  of  the  solved  equation  is 

dc  sc  sc  d  (n  ac)  a  (n  ac\ 

dt  dx  dy  dx  [  x  dx)  dy  {  y  dy )  (a5) 

+  aj  C2  +  a  =  0 


where 

C  =  concentration  of  sediment 
u  -  depth-integrated  velocity  in  x-direction 
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v  =  depth-integrated  velocity  in  y-direetion 
Dx  -  dispersion  coefficient  in  x-direction 
Dy  -  dispersion  coefficient  in  y-direction 
a.j  =  coefficient  of  concentration-dependent  source/sink  term 
a2  =  coefficient  of  source/sink  term 

The  source/sink  terms  in  Equation  B5  are  computed  in  routines  that  treat 
the  interaction  of  the  flow  and  the  bed.  Separate  sections  of  the  code  handle 
computations  for  clay  bed  and  sand  bed  problems. 


Sand  transport 

The  source/sink  terms  are  evaluated  by  first  computing  a  potential  sand 
transport  capacity  for  the  specified  flow  conditions,  comparing  that  capacity 
with  amount  of  sand  actually  being  transported,  and  then  eroding  from  or 
depositing  to  the  bed  at  a  rate  that  would  approach  the  equilibrium  value  after 
sufficient  elapsed  time. 

The  potential  sand  transport  capacity  in  the  model  is  computed  by  the 
method  of  Ackers  and  White  (1973),  which  uses  a  transport  power  (work  rate) 
approach.  It  has  been  shown  to  provide  superior  results  for  transport  under 
steady-flow  conditions  (White,  Milli,  and  Crabbe  1975)  and  for  combined 
waves  and  currents  (Swart  1976).  Flume  tests  at  the  U.S.  Army  Engineer 
Waterways  Experiment  Station  have  shown  that  the  concept  is  valid  for 
transport  by  estuarine  currents. 

The  total  load  transport  function  of  Ackers  and  White  is  based  upon  a 
dimensionless  grain  size 


g(s  ~  1) 


1/3 


(A6) 


where 

D  =  sediment  particle  diameter 
s  =  specific  gravity  of  the  sediment 
v  =  kinematic  viscosity  of  the  fluid 

and  a  sediment  mobility  parameter 
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where 

x  =  total  boundary  shear  stress 

n'  =  a  coefficient  expressing  the  relative  importance  of  bed-load  and 
suspended-load  transport,  given  in  Equation  A9 
t’  =  boundary  surface  shear  stress 

The  surface  shear  stress  is  that  part  of  the  total  shear  stress  which  is  due  to  the 
rough  surface  of  the  bed  only,  i.e.,  not  including  that  part  due  to  bed  forms 
and  geometry.  It  therefore  corresponds  to  that  shear  stress  that  the  flow  would 
exert  on  a  plane  bed. 

The  total  sediment  transport  is  expressed  as  an  effective  concentration 
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P  A  h  [n1  t 

(A8) 

where  V  is  the  average  flow  speed,  and  for  1  <  D„  <  60 

o' 

ri  =  1.00  -  0.56  log  Dgr 

(A9) 

A-  °-23  +  0.14 

K 

(A10) 

log  C  =  2.86  log  Dgr  -  (log  Dgr)2  -  3.53 

(All) 

it.  -  966  ♦  1.34 

(A  12) 

For  <  60 

gf 

n’  =  0.00 

(A13) 

A  =  0.17 

(A14) 
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C  «  0.025 


(A15) 


m  =  1.5 

(A16) 


Equations  A6-A16  result  in  a  potential  sediment  concentration  Gp.  This 
value  is  the  depth-averaged  concentration  of  sediment  that  will  occur  if  an 
equilibrium  transport  rate  is  reached  with  a  nonlimited  supply  of  sediment 
The  rate  of  sediment  deposition  (or  erosion)  is  then  computed  as 


(A17) 


where 

C  =  present  sediment  concentration 
tc  =  time  constant 

For  deposition,  the  lime  constant  is 


tc  -  larger  of 


At 

or 

Cdh 
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and  for  erosion  it  is 


tc  =  larger  of  j 


At 

or 

CJ 

U 


(A  19) 


where 

At  =  computational  time-step 

Crf  -  response  time  coefficient  for  deposition 

Vs  =  sediment  settling  velocity 

Ce  =  response  time  coefficient  for  erosion 

The  sand  bed  has  a  specified  initial  thickness  which  limits  the  amount  of 
erosion  to  that  thickness. 
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Cohesive  sediments  transport 


Cohesive  sediments  (usually  clays  and  some  silts)  are  considered  to  be 
depositional  if  the  bed  shear  stress  exerted  by  the  flow  is  less  than  a  critical 
value  x When  that  value  occurs,  the  deposition  rate  is  given  by  Krone’s 
(1962)  equation 
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where 

S  =  source  term 

=  fall  velocity  of  a  sediment  particle 
li  =  flow  depth 

C  =  sediment  concentration  in  water  column 
x  =  bed  shear  stress 
xd  =  critical  shear  stress  for  deposition 
Cc  -  critical  concentration  =  300  mg/V 

If  the  bed  shear  stress  is  greater  than  the  critical  value  for  particle  erosion 
xe,  material  is  removed  from  the  bed.  The  source  term  is  then  computed  by 
Ariathurai’s  (Ariathurai,  MacArthur,  and  Krone  1977)  adaptation  of 
Partheniades’  (1962)  findings: 
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where  P  is  the  erosion  rate  constant,  unless  the  shear  stress  is  also  greater  than 
the  critical  value  for  mass  erosion.  When  this  value  is  exceeded,  mass  failure 
of  a  sediment  layer  occurs  and 
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where 

Tl  -  thickness  of  the  failed  layer 
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All 


PL  =  density  of  the  failed  layer 
At  =  time  interval  over  which  failure  occurs 
=  bulk  shear  strength  of  the  layer 

The  cohesive  sediment  bed  consists  of  1  to  10  layers,  each  with  a  distinct 
density  and  erosion  resistance.  The  layers  consolidate  with  overburden  and 
time. 


Bed  shear  stress 

Bed  shear  stresses  are  calculated  from  the  flow  speed  according  to  one  of 
four  optional  equations:  the  smooth-wall  log  velocity  profile  or  Manning 
equation  for  flows  alone;  and  a  smooth  bed  or  rippled  bed  equation  for  com¬ 
bined  currents  and  wind  waves.  Shear  stresses  are  calculated  using  the  shear 
velocity  concept  where 


xb  =  P“? 
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where 

Tb  =  bed  shear  stress 
u%  =  shear  velocity 

and  the  shear  velocity  is  calculated  by  one  of  four  methods: 
a.  Smooth-wall  log  velocity  profiles 
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which  is  applicable  to  the  lower  15  percent  of  the  boundary  layer  when 


v 


where  u  is  the  mean  flow  velocity  (resultant  of  u  and  v  components) 
b.  The  Manning  shear  stress  equation 
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where  CME  is  a  coefficient  of  1  for  SI  (metric)  units  and  1 .486  for 
non-SI  units  of  measurement. 

c.  A  Jonsson-type  equation  for  surface  shear  stress  (plane  beds)  caused  by 
waves  and  currents 
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where 

fw  =  shear  stress  coefficient  for  waves 
uom  ~  max‘mum  orbital  velocity  of  waves 
fc  =  shear  stress  coefficient  for  currents 

d.  A  Bijker-type  equation  for  total  shear  stress  caused  by  waves  and 
current 
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Solution  method 

Equation  A5  is  solved  by  the  finite  element  method  using  Galerkin 
weighted  residuals.  Like  RMA-2V,  which  uses  the  same  general  solution  tech¬ 
nique,  elements  are  quadrilateral  and  may  have  parabolic  sides.  Shape  func¬ 
tions  are  quadratic.  Integration  in  space  is  Gaussian.  Time-stepping  is  per¬ 
formed  by  a  Crank-Nicholson  approach  with  a  weighting  factor  (0)  of  0.66.  A 
front-type  solver  similar  to  that  in  RMA-2V  is  used  to  solve  the  simultaneous 
equations. 
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